Experimental and numerical study of steady state stability in a toluene biodegrading biofilter

Different steady states in a toluene biodegrading biofilter were explored experimentally and numerically. Experimental results showed that a gradual increase of the toluene inlet concentration over several weeks leads to a consistently low exit concentration, with a drastic increase at an inlet concentration change from 7.7 to 8.5 g m−3, indicating an alteration in steady state. A significant and sudden drop in the removal efficiency from 88 to 46% was observed. A model that includes nitrogen and biomass dynamics predicted results matching the experimental biofilter performance well, but the timing of the concentration jump was not reproduced exactly. A model that assumes a gradual increase of toluene inlet concentration of 0.272 g m−3 per day, accurately reproduced the experimental relationship between inlet and outlet concentration. Although there was variation between experimental and simulated results, a clear confirmation of the jump from one steady state to another was found.

Gaseous emissions from various industrial processes need to be treated before they are released into the environment because of their harmful effect on humans, animals, and the environment. Waste gas streams potentially contain hazardous substances such as volatile organic compounds (VOCs) or odorous compounds. The oil and gas, chemical and pharmaceutical industries are anthropogenic sources of such VOCs. The need to comply with stringent legislation can pose a challenge, because any waste gas treatment system must be reliable, sustainable, and efficient. Physical and chemical methods are applied in industry but are often associated with disadvantages such as operational costs as well as the generation of undesirable by-products compared to their biological counterpart [1][2][3][4] . The treatment of polluted air by means of biological techniques is known and well-established [5][6][7][8][9][10][11] . In a biofilter, microbes degrade the compound of interest and produce water, carbon dioxide and heat 2,4 . A biofilter simulation enables the prediction of the behavior and efficiency of a system and is a substantial part of biological air pollution control research in recent decades. It is an auxiliary method to aid biofilter designers and manufacturers, as they link biological growth kinetics to reaction engineering.
Ottengraf and Van Den Oever 12 developed one of the first models in 1983. In their study, experimental results corresponded well with a model that assumes mass transfer limitation and zero-order kinetics. Since then, numerical based approaches, sometimes supported by experimental data of biofilter and biofilm kinetics were developed in the last decades [13][14][15][16][17][18][19][20] . These models overcome some of the limitations of the Ottengraf and Van den Oever model. Yan et al. 21 studied the effect of the filter bed porosity on removal efficiency. They showed a low impact (i.e., high efficiency) at low Darcy numbers (< 10 -4 ), whereas at high Darcy numbers, the porosity significantly affects the removal efficiency. Woudberg et al. 22 proposed a pore-scale model to predict the pressure drop of a biofilter. The prediction of Malakar et al. 23 returned to the original Ottengraf and Van Den Oever model to describe the theoretical elimination capacity and theoretical average biofilm thickness in a toluene biofilter. The biofilm was assumed to be static in this study. Dorado et al. 24 estimated the kinetic parameters of a more sophisticated model involving biofilm diffusion and Monod kinetics. Biofilm thickness was also assumed to be static in this study. With continuous digitalization and increasing available biofiltration data, the implementation of machine learning approaches to predict the biofilter performance is now possible 25 .
In Süß and De Visscher 26 , it was shown that a sudden drop in removal efficiency (RE) potentially is an implication of different steady states occurring in the biofilm if (1) the pollutant concentration exceeds a threshold value and (2) substrate degradation and substrate inhibition follow Haldane kinetics. However, the dynamics and implication of steady states in a biofilm or biofilter aerobically treating gaseous contaminants are not yet fully understood, mainly because existing biofilter models do not include the dynamic nature of biofilm growth and die-off during biofiltration. In this study, we hypothesize that biofilters can progress from a single steady state behavior to a multi steady state behavior, driven by the growth of the biofilm thickness as the biofilter is exposed to increasing pollutant concentrations. At the early stare of biofilter operation, biofilms are too thin to display multiple steady states. As biofilms grow thicker, the filter goes through a phase of high activity at conditions where a second, low-activity steady state exists. At this stage, a slight concentration increase can cause the biofilter to "crash", i.e., to drop from the high activity to the low activity steady state.
To test this hypothesis, a toluene biofiltration experiment was run with inlet concentrations gradually increasing to high levels (> 10 g/m 3 ), and a model was developed with a dynamic biofilm thickness driven by nitrogen availability for biomass to grow on toluene.

Materials and methods
Microorganism and filter bed material of biofilter. A sterilized mixture of compost and wood chips (vol% 80/20) was used as filter bed. In Table 1, the filter bed analysis is shown. Analysis was conducted by AGAT Laboratories. Chemicals such as toluene and solids for the cultivation medium were obtained from Sigma Aldrich. An adapted liquid BH-medium 17 was used as a growth medium and contained 1 g L −1 KH 2 PO 4 , 1 g L −1 www.nature.com/scientificreports/ Na 2 HPO 4 , 0.5 g L −1 NH 4 NO 3 , 0.002 g L −1 FeCl 3 , 0.002 g L −1 MnSO 4 ·2H 2 O, 0.2 g L −1 MgSO 4 ·7H 2 O and 0.02 g L −1 CaCl·2H 2 O. The medium was sterilized before use. Toluene was used as a single carbon source. The growth medium was used to culture a toluene degrader, Nocardia sp. 17 . To introduce additional microorganisms into the inoculant, an air stream was piped into the serum bottles. After 3 days, plating of cells was conducted, and standard procedures were followed 27 . Two other bacterial strains were found. These other two bacteria were tested on their ability to biodegrade toluene as a single carbon source by firstly isolating and secondly growing them in serum bottles with Toluene as only carbon source. No microbial growth was determined hence, these microbes were not able to biodegrade toluene. Therefore, it is assumed that Nocardia sp. is the only Toluene-degrading microbe in the mixture and that the other strains were autotrophic bacteria. No further characterization of microbes was conducted. Nocardia sp. as Toluene-degrading and two non-Toluene degrading microbes were cultivated under ambient temperature (21-22 °C) in serum bottles sealed with a butyl rubber septum and aluminum crimp.
Biofilter set-up and experimental conditions. The Toluene biofiltration experiments were carried out in a lab-scale biofilter with a packing as described above, inoculated with 100 mL of inoculant (OD 650 around 0.2). The column was made of polycarbonate with a total height of 63 cm and an internal diameter of 10 cm. At the bottom, glass beads at a height of 5 cm and a 1 cm perforated Plexiglas plate were used to evenly distribute the inlet gas stream. Sampling ports to measure gas samples were centered at the top and bottom of the column and sealed with GC septa (0.95 cm diameter). The biofilter was filled with sterilized filter media, containing dispersed inoculant to a height of 25 cm and placed in a fume hood under ambient temperatures of 21-22 °C. The experimental set-up is based on 17 . An air pump (pond master Ap-40) was used to generate one air stream, which was split and sent through two gas washing bottles, one filled with tap water (resistivity 0.0029 MΩ, hardness 169 mg/L as CaCO 3 , pH 7.8), and one filled with the Toluene. Subsequently, the two gas flows were combined and mixed in an empty gas-washing bottle before being introduced into the biofilter. The flow was controlled by two rotameters (Cole-Parmer) located after the gas flow split, and total flow rate was measured (TI-400) before the sample port and inlet. Biofilter operation. An air stream contaminated with toluene was treated under different operational conditions. An EBRT of 4.5 min was maintained and a stepwise increase of the inlet concentration was conducted. This EBRT is higher than typical but not unusual when the inlet concentration is high [28][29][30] . An adsorption test for toluene was conducted with sterilized filter packing and indicated no adsorption.
The operational parameters of the biofilter experiments are shown in Table 2.
Analytical methods. An SGE 250 µL gastight syringe was used to draw and inject 200 µL of gas sample from the gas sampling ports and into the analyzer. The gas samples were analyzed with a gas chromatograph (GC-2014, Shimadzu) equipped with an FID and Rtx®-Wax capillary column (30 m × 0.53 mm × 1 µm). The injector and detector temperatures were set at 250 °C. The oven temperature for toluene was 80 °C. Helium was used as a carrier gas.
Model description. The multiplicity of steady states in an aerobic biofilter is rarely studied and understood.
Simulations with a biofilm model by Süß and De Visscher 26 showed that different steady states can exist inside a biofilm, considering no other back-mixing mechanisms. It was also established that a change in steady state can lead to a sudden decrease in RE for a small increase of the pollutant concentration. To further investigate the possible steady states in a biofilter, a full biofilter model with biofilm dynamics was developed based on a number of assumptions about the biofilter and the biofilm: 1. Gas phase flow is assumed to behave in plug flow pattern. Thus, axial dispersion is neglected.
2. Due to low gas phase concentration, the gas-biofilm equilibrium at the interface is described by Henry's law. Model development. To minimize computation time of the model, orthogonal collocation 31 was used both in the gas phase and in the biofilm phase, linked together by using Henry's law at the biofilm surface. In order to predict the outlet concentration of the biofilter, the packing material was divided into 25 collocation points along the biofilter height. For each collocation point in the gas phase, the gas-phase concentration was calculated, and a biofilm consisting of 10 collocation points was modeled to calculate the average reaction rate, net growth rate, and biofilm concentration profile. In addition, a nitrogen cycle was also considered in the model 32  In order to calculate the concentration profile in the biofilm, diffusion and reaction rate, Equations (5) and (6), were linked together, considering the biofilm thickness L [m], the distance coordinate in the biofilm x [m] and a dimensionless distance coordinate in the biofilm x′ (= x/L) using a material balance, which leads to the following expression. where the inside boundary of the biofilm away from the gas is represented by x = 0. The partial differential equation in Eq. (7), which is first-order in time, and second-order in space, was solved by Orthogonal Collocation to approximate the concentration profile in the biofilm and as outlined by Villadsen and Stewart 31 . Next, the concentration in the gas phase in each collocation point was computed by equating the transfer of toluene in the gas phase towards the biofilm to the integrated toluene biodegradation in the biofilm. To account for the biofilm growth, the net growth rate of the microorganism needs to be considered. As mentioned above, the nitrogen cycle is considered in this model and a part of it is expressed in the following equation: where, µ max , µ net , N inorg k N_Nitrogen  where f N is the nitrogen fraction of the microorganisms, and µ is the average growth rate over the biofilm, calculated in a manner similar to r . A is the specific area of the biofilm [m 2 kg compost −1 ]. The dynamics of inorganic nitrogen is calculated as follows and shown in 32 : where N inorg , k minN and k uptakeN represent the inorganic nitrogen content in the packing material [g N kg dw −1 ], the nitrogen mineralization rate constant [h −1 ], and the nitrogen uptake rate constant [h −1 ], respectively.
The thickness of the biofilm, L, is calculated with the following equation:

Results and discussion
Experimental data-steady state. In Fig. 1 the inlet and corresponding outlet concentrations of the main experiment, are displayed to evaluate the occurrence of two steady states. The standard deviation of the exit concentration measurements was 0.10 g m −3 (i.e., a standard error of 0.056 g m −3 based on triplicate measurements). Until day 31, the outlet concentration barely changed, although the inlet concentration was gradually increased. At this time a decline of the RE from 99 to 88% was measured as the inlet concentration was increased by 4.9 g m −3 (increase from 2.9 to 7.7 g m −3 ). The inlet concentration was adjusted and increased prior to day Scientific Reports | (2022) 12:12510 | https://doi.org/10.1038/s41598-022-15620-w www.nature.com/scientificreports/ 31 to allow the biomass to adjust to the new load prior to the measurement. This is indicated in Fig. 1, 3 and 4 with a theoretical calculated inlet concentration depicted as a rhomboid. The next increase of the inlet concentration was by 0.8 g m −3 (to a value of 8.5 g m −3 ), which led to a steep increase of the outlet concentration. This corresponds to a decrease of RE from 88 to 46%. Such a significant decline in RE for a small increase in inlet concentration could indicate that a boundary of a stable steady state region has been crossed. As discussed and numerically shown in Süß and De Visscher 26 with a biofilm model, a change in steady state can be explained by substrate degradation and substrate inhibition following Haldane kinetics and the diffusion behavior of the pollutant into the biofilm. Distinctive for Haldane kinetics is that low reaction rates occur at low and high concentrations and high reaction rates occur at medium concentrations. With regard to the conducted experiments, the transition from high to low reaction rates is in the range of 7.7-8.5 g m −3 . Furthermore, diffusion limitations are an important factor as well. When medium range concentrations are present at the surface of the biofilm, it is possible to maintain such concentrations throughout the biofilm, leading to a high reactivity and pronounced diffusion limitation. Hence, a significant concentration gradient will be upheld in the biofilm and will consequently maintain a medium range concentration at the inside of the biofilm. On the other hand, if a high enough concentration is present at the surface of the biofilm, reaction rates near the surface are low. Thus, a high concentration can develop throughout the biofilm and therefore result in low reactivity. In this case, diffusion limitation is not pronounced.
Fitting computer simulation to experimental data-steady states. The computer simulation developed here was used to predict the outlet concentrations based on obtained inlet concentrations of the experimental trial. To optimize the model, parameters A and k min were used as adjustable variables. The remaining parameter values were taken from previous studies 33,34 . Used model parameters are displayed in Table 2 and the simulation results and experimental data are shown in Fig. 2.
The model prediction and the experimental results are in good agreement until day 23. After that day, the predicted outlet concentration suddenly rose and declined over a period of 8 days. As indicated above, the inlet concentration was increased on day 24, and the exit concentration measurement was made on day 31. This is indicated in Fig. 2, with the expected inlet concentration after adjustment depicted as a rhomboid. That increase of the inlet concentration caused the predicted sudden rise at the outlet. The subsequent predicted decline of outlet concentration is possible due to the adaption of the system to the high stepwise increase of the inlet concentration (biofilm growth). During the next time period (day 31 to 37) the experimental trial indicates a change in steady states, whereas the model prediction indicates such a change between day 37 and 38, at a concentration change from 8.502 to 9.257 g m −3 . This corresponds to a 47% decline in RE. Between day 37 and 40 the maximum achieved decline in RE is 72%. The root mean square error (RMSE) between the modeled and measured outlet toluene concentration is 1.39 g m −3 (r 2 = 0.900). Although the change in steady state is not predicted in the same time period and concentration range, an indication of change can be seen. A further adjustment of model parameters could possibly increase the accuracy of the model. Alternatively, it may be the case www.nature.com/scientificreports/ that the nitrogen dynamics model of Eqs. (15)(16) does not fully capture the dynamics of biofilm development, particularly under rapidly changing conditions. If this is the case, then it can be expected that a model with a more gradually increasing inlet concentration can describe the relationship between inlet and outlet toluene concentration more satisfactorily.
To test this hypothesis, a second model run was carried out, where the inlet concentration was steadily increased with uniform increments of 0.272 g m −3 per day. The results are shown in Fig. 3, overlaid with the experimental data. To optimize the model fit, the surface area (A) was increased slightly from 0.95 to 1.1 m 2 kg −1 .   Fig. 3 on day 30 for the model and day 37 for the experiment. Consequently, in both cases, simulation and experimental results, a change in steady state occur at similar concentration. The modeled jump is sharper than the observed jump. This is because the model assumes cross-sectionally uniform biofilm thickness, whereas the actual biofilm will not be uniform within a cross-section. The inlet concentration before the jump was 7.909 g m −3 for the simulation and 7.7 g m −3 for the experiment, and the inlet concentration after the observed jump were 8.181 g m −3 and 8.5 g m −3 for the simulation and experiment, respectively. When expressed as a function of time, the RMSE between the model and the experimental results is 1.36 g m −3 (r 2 = 0.907). Figure 4 is included to better illustrate the similarity between the experimental data and the simulation shown in Fig. 3. In Fig. 4 the outlet concentration is plotted versus inlet concentration. A very good agreement is obtained (RMSE = 0.40 g m −3 , r 2 = 0.992). In this figure, a significant increase in the outlet concentration is observed at a small change in inlet concentration.
This indicates a jump from a high activity steady state to a low activity steady state. The observed behavior should not be confused with a collapse of the biofilm. A collapsing biofilm at such a low inlet concentration change would occur over an extended time (days to weeks), not suddenly, considering no other changing factor or inhibition.
In the above simulation, the concentration in the biofilm is low when high activity is predicted and high when low activity is shown. This indicates a non-saturated and saturated biofilm, respectively. This is illustrated in Fig. 5, where the concentration profile in the biofilm is shown at the beginning of day 37 (888 h) and one hour later (889 h), in biofilter grid point 6 (i.e., at a biofilter height of 3.25 cm), simulating the experimental conditions. The concentration is plotted as a function of distance from the solid surface of the packing material to the surface of the biofilm. As can be seen, the concentration sharply declines towards the inside of the biofilm on day 37, which indicates a non-saturated biofilm and high activity. One hour later, on day 37.04, the concentration in the biofilm has settled to a new steady state associated with a nearly saturated biofilm and lower activity.
Fitting computer simulation to experimental data-biofilter. To further validate the accuracy and applicability of the computer simulation, it was verified against a second set of biofilter experiments 33 . In this experiment, the inlet concentration was increased stepwise and held at each stage until a steady state was reached. In Fig. 6 the predicted outlet and experimental inlet and outlet concentrations are shown, and the model parameters are listed in Table 2. The experimental outlet concentrations have a standard deviation of 0.071 g m −3 (standard error of triplicates 0.041 g m −3 ). To reach a good fit between the model and the data, the Outlet concentration [g m -3 ] Inlet concentration [g m -3 ] expeimental data simulaƟon www.nature.com/scientificreports/ specific surface area was reduced by a factor 4. This corresponds with increasing the packing size from about 6 mm to about 25 mm (assuming spherical particles and a solid density of 1000 kg m −3 ). These are reasonable values, and the increase was expected because the straw used as a bulking agent in these experiments did not sustain the structure of the biofilter as well as the wood chips used in the first experiment. The model follows the experimental data well until day 36. After that, the model underpredicts the outlet concentration. It is hypothesized that settling of the biofilter reduced the biofilm specific surface area exposed Concentration [g m -3 ]

Time [days]
Inlet experiments Outlet experiments simulaƟon www.nature.com/scientificreports/ to the gas phase. Some slight settling was observed in this biofilter, unlike the biofilter of the data in Fig. 1. The reduced specific surface area reduced the biofilm area in the biofilter. Despite the lack of fit towards the end of the experiment, the model showed a RMSE of only 0.153 g m −3 (r 2 = 0.937). When only the first 36 days are considered, the RMSE is 0.059 g m −3 (r 2 = 0.950). The simulations indicate that the biofilm is strongly diffusionlimited in this biofilter, so a direct proportionality between biofilm area and activity is expected. The assumption of proportionality between activity to specific surface area is consistent with Delhoménie et al. 35 , who found that maximum EC decreased with increasing particle size, but increases with increasing specific surface area. This implies that at that point in time the biofilm surface area has a more pronounced impact on the RE then the activity of the biofilm. These findings are closely linked to the moderately hydrophobic nature of toluene. Zhu et al. 36 found that a toluene biofilter has an degradation efficiency intermediate between isobutanol (H = 0.0005) and n-hexane (H = 53), indicating partial diffusion limitation, consistent with our model interpretation of the results. Nevertheless, biofiltration can be effective at Henry's law constants of 10 and above (Haque et al., 37 ). Ranjbar and Gheeni 38 found that the Henry constant and the specific surface area were the most sensitive parameters in their biofilter model, consistent with our findings. Kalantar et al. 39 found that toluene is more strongly diffusion-limited in a two-phase biotrickling filter than was found in our work. However, biotrickling filters have thicker biofilms than biofilters, which explains the difference.
The accuracies of the model predictions presented here are similar to the accuracies that are typical for biofilter models, such as Ranjbar and Ghaemi 38 , and San Valeo et al. 40,41 .

Conclusion
Experimental and simulation results of a toluene biofilter with a steadily increasing inlet concentration show a jump from a high-to low-activity steady state, albeit at a slightly different timing. Results showed a good overall prediction of the outlet concentration, except at the end of the experiment, where settling of the filter bed material may have reduced the biofilm area. An investigation of modeled toluene concentration profile in the biofilm before and after the sudden jump in RE confirmed that the cause of the jump is a transition from a diffusionlimited high-activity state to a low-activity state.

Data availability
The datasets used and/or analysed during the current study available from the corresponding author on reasonable request.